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An extension of the Asakura-Oosawa-Vrij model of hard sphere colloids and non-adsorbing poly- 
mers, that takes polymer non-ideality into account through a repulsive stepfunction pair potential 
between polymers, is studied with grand canonical Monte Carlo simulations and density functional 
theory. Simulation results validate previous theoretical findings for the shift of the bulk fiuid demix- 
ing binodal upon increasing strength of polymer-polymer repulsion, promoting the tendency to mix. 
For increasing strength of the polymer-polymer repulsion, simulation and theory consistently pre- 
dict the interfacial tension of the free colloidal liquid-gas interface to decrease significantly for fixed 
colloid density difference in the coexisting phases, and to increase for fixed polymer reservoir packing 
fraction. 
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I. INTRODUCTION 



Various different levels of description have been em- 
ployed in order to study mixtures of colloidal parti- 
cles and non-adsorbing globular polymers P, Q • While 
colloid-colloid interactions are reasonably well described 
by those of hard spheres, the effective interactions be- 
tween colloids and polymers, as well as that between the 
centers of two polymers, are more delicate. The model 
due to Asakura and Oosawa |3] and Vrij 4] (AOV), tak- 
ing the colloids to be hard spheres and the polymer- 
polymer interactions to vanish, is arguably the most sim- 
ple description of real colloid-polymer mixtures. Its ap- 
peal clearly lies in its simplicity, rather than in very close 
resemblance of experimental colloid-polymer systems, en- 
abling detailed simulation (sj, |^J3JE_Pl and theoretical 

trapjlljllstudies of bulk ifllMlg and interfacial 
H EHUHEB 111 m properties. More realistic ef- 
fective interactions between the constituent particles can 
be systematically obtained by starting at the polymer 
segment level, and integrating out the microscopic de- 
grees of freedom of the polymers 's', . The resulting 
colloid-polymer interaction is a smoothly varying func- 
tion of distance, and excluded volume between polymer 
segments leads to a soft Gaussian-like polymer-polymer 
repulsion. Such effective interactions have been used to 
calculate bulk phase behavior |^l2l|. 

In order to retain most of the simplicity of the AOV 
model, but still to capture polymer non-ideality, the AOV 
model was extended using a repulsive stepfunction pair 
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potential between polymers in Ref. |23. The colloid- 
polymer interaction was kept as that of hard spheres. 
This introduces one additional model parameter, namely 
the ratio of stepheight and thermal energy, interpolating 
between the AOV case when this quantity vanishes, and 
the additive binary hard sphere mixture when it becomes 
infinite. Following the hard sphere and AOV cases 
[r2i. ..13j . the well-defined particle shapes of this model 
allowed to obtain a geometry-based density functional 
theory (DFT) [Mill specifically tailored for this model 
22] . The theory can, in principle, be applied to arbitrary 
inhomogeneous situations, which constitutes an a poste- 
riori justification for using these model interactions. The 
trends found for the fluid-fluid demixing transition into a 
colloidal liquid phase (that is rich in colloids and poor in 
polymers) and a colloidal gas phase (that is poor in col- 
loids and rich in polymers), upon increasing the polymer- 
polymer interaction strength, demonstrated improved ac- 
cordance with the experimental findings of Ref. .26 . as 
compared to the AOV case 0, , where the DFT pre- 
dicts the same phase diagram as the free volume theory 
Furthermore the extended AOV model was used in 
the study of the "floating liquid phase" that was found in 
sedimentation-diffusion equilibrium p^. Further discus- 
sion of the model and its relation to the AOV prescription 
is given in Ref. |2^ 

The aim of the present paper is twofold. First, we want 
to assess in detail the accuracy of the DFT of Ref. [23 by 
comparing to results from direct simulations of the ex- 
tended AOV model. While phase-coexistence is often 
studied in the Gibbs ensemble [S^, we instead choose 
to take advantage of recently developed methods 0, 
that rely on the grand canonical ensemble. Benefits are 
accurate estimates of the interfacial tension 7, 30] and 
access to the critical region |3l|. The binodal we ob- 
tain in the simulations is compared to that from DFT, 
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and good agreement between both is found. This result 
strongly supports the original claim f23| that a straight- 
forward perturbation theory, taking the AOV system as 
the reference state and adding polymer-polymer interac- 
tions in a mean-field like manner, will fail, as the bulk 
fluid-fluid binodal from this approach differs markedly 
from that of the geometry-based DFT, as discussed in de- 
tail in Ref. 22. Instead, even to lowest order in strength of 
the polymer-polymer interaction, a non-trivial coupling 
to the colloid density field needs to be taken into ac- 
count; the DFT of Ref. |22 does this intrinsically. Our 
second aim is to study the interfacial tension between 
demixed colloidal gas and liquid phases, which is known 
from experiments pi El 13 El , theory Hi El [13, El 
and simulation 0, l37| | to be orders of magnitude smaller 
than that of atomic substances. While much work has 
been devoted to the case of non-interacting polymers 
|3lillIllI3llllIl|3l|33, only quite recent studies 
addressed the effect of polymer non-ideality 38, 39, 4^ . 
Aarts et al. use a square gradient approach based 
on the free energy of a free volume theory, to include 
polymer interactions j4H . and the mean spherical ap- 
proximation for the direct correlation function. They find 
the gas-liquid interfacial tension to decrease as compared 
to the case of non-interacting polymers, when plotted 
as a function of the density difference in the coexisting 
phases. Similar findings were reported by Moncho-Jorda 
et al. 39] , who also used a square gradient approach, but 
based on an effective colloid-colloid depletion potential 
that reproduces simulation results accurately (42]. Our 
present study goes beyond Refs. IH IH and also beyond 
the very recent Ref. l40l as we employ a non-perturbative 
DFT treating the full two-component mixture of colloids 
and interacting polymers. We compare results for the in- 
terfacial tension to our data from direct simulation of the 
binary mixture. Besides its intrinsic interest, the inter- 
facial tension is moreover relevant for the occurrence of 
capillary condensation in confined systems 43, 44, 45] as 
is apparent from a treatment based on the Kelvin equa- 
tion [i^. 

The paper is organized as follows. In Sec. ^ we define 
the extended AOV model taking into account polymer- 
polymer repulsion. In Sec. IIIII we briefly sketch the the- 
oretical and simulation techniques. In Sec. lIVI results for 
fiuid-fiuid phase behavior and the interfacial tension are 
presented, and we conclude in Sec. El 



center distance r between two particles, given as 




oo r < (Jc 

otherwise, 

oo r < (dc + fp)/2 

otherwise, 

e r < (Tp 

otherwise. 



(1) 
(2) 
(3) 



Particle numbers are denoted by iV^, and as bulk ther- 
modynamic parameters we use the packing fractions ?/; = 
TTcrfNi/ {6V), where V is the system volume. As an alter- 
native to 77p, we use the packing fraction rj^ — 7r(TpPp/6 
in a reservoir of pure polymers, interacting via V^pp(?') as 
given above, that is in chemical equilibrium with the sys- 
tem, with Pp being the polymer number density in the 
reservoir. 

The model is characterized by two dimensionless con- 
trol parameters, namely the polymer-to-colloid size ratio 
q = CTp/cTc, and the scaled strength of polymer-polymer 
repulsion f3e, with (3 — l/{kQT), kn the Boltzmann con- 
stant, and T the absolute temperature. As limiting cases, 
the present model possesses the AOV model for /3e = 0, 
where polymer-polymer interactions are ideal, and the bi- 
nary additive hard sphere mixture in the limit /?e — > oo. 
One can use the parameter (3e to match to a real system 
at a given thermodynamic statepoint, polymer type and 
solvent, by imposing equality of the second (polymer- 
polymer) virial coefficients. See Ref. for an in-depth 
discussion of this procedure. 



III. METHODS 

A. Density functional theory 

To investigate bulk and interfacial properties of the 
present model, we use the geometrically-based DFT of 
Ref. I23 that has its roots in generalizations of Rosenfeld's 
fundamental-measure theory for additive hard sphere 
mixtures |23j ]. namely the treatment of the "penetrable 
sphere model" ^3] (equivalent to the present polymer 
particles) and the DFT genuinely developed for the AOV 
model |l2l ]. The minimization of the grand potential is 
carried out with a simple iteration technique. We refer 
the reader directly to Ref. 22i for all details about the 
(approximate) Helmholtz free energy functional. 



II. MODEL 



Simulation method 



We consider a mixture of hard sphere colloids (species 
c) of diameter CTc, and effective polymer spheres (species 
p) of diameter cTp, that interact via pairwise potentials 
Vij{r) with (i,j) £ (c, p), as function of the center-to- 



The simulations are carried out in the grand canoni- 
cal ensemble, where the fugacities Zc and Zp, of colloids 
and polymers, respectively, and the total volume V are 
fixed, while the numbers of particles, Nc and A^p, are 
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allowed to fluctuate. We use a rectangular box of dimen- 
sions L X L X D, with periodic boundary conditions in all 
three directions, and simulate the full mixture as defined 
by the pair potentials given in Eqs. Q-®, i-C- the posi- 
tional degrees of freedom of both colloids and polymers 
are explicitly taken into account. Note that asymmetric 
binary mixtures are in general difficult to simulate, and 
prone to long equilibration times. To alleviate thisprob- 
lem, we rely on a recently developed cluster move 0, 13 ' 
that has already been applied successfully to the stan- 
dard AOV model 0, ISlL I4M l49j| . Here we perform the 
generalization to /3e > which is straightforward. 

During the simulation, we measure the probability 
P(?7c) of finding a certain colloid packing fraction ijc- 
At phase coexistence, the distribution P{r]c) becomes bi- 
modal, with two peaks of equal area, one located at small 
values of 77c corresponding to the colloidal gas phase, 
and one located at high values of 7]c corresponding to 
the colloidal liquid phase. Typical coexistence distri- 
butions for the standard AOV model can be found in 
Ref.Q, and our present results display similar behavior. 

The equal area rule l50l| implies that /q''''^ P{'i]c)'iVc = 
^ P{ric)dric, with (rjc) the average of the full distri- 
bution (rjc) = ilcP{Vc)'i''lc, where we assume that 
P{r]c) has been normalized to unity P{r]c)dric — 1. 
The packing fraction of the colloidal gas rj^ now fol- 
lows trivially from the average of P{tIc) in first peak 
?7f = 2 /j''"^ r]cP{r]c)dr]c, and similarly for the colloidal 
liquid ??' = 2 ^ r]cP{ric)dr]c, where the factors of 2 are 
a consequence of the normalization of P{r]c)- 

The interfacial tension 7 is extracted from the loga- 
rithm of the probability distribution W = k^T In P{r]c). 
Note that —W corresponds to the free energy of the sys- 
tem. Therefore, the height of the peaks in W, mea- 
sured with respect to the minimum in between the peaks, 
equals the free energy barrier separating the colloidal gas 
from the colloidal liquid "s^. Fl may be related to the 
interfacial tension 7 by noting that, at colloid packing 
fractions between the peaks ryf -C 77c vl^ ttio system 
consists of a colloidal liquid in coexistence with its va- 
por. A snapshot of the system in this regime, would 
reveal a so-called slab geometry, with one region dense 
in colloids, and one region poor in colloids, separated by 
an interface (because of periodic boundary conditions, 
two such interfaces are actually present). If an elongated 
simulation box with D > L is used, rather than a cubic 
box with D ^ L, the interfaces will be oriented perpen- 
dicular to the elongated direction, since this minimizes 
the interfacial area, and hence the free energy of the sys- 
tem. The total interfacial area in the system thus equals 
2L2 and, following Ref.lM 7 = FlI{2L'^). An additional 
advantage of using an elongated simulation box is that 
interactions between the interfaces are suppressed. This 
enhances a flat region in W between the peaks, which 
is required for an accurate estimate of the interfacial 
tension. In this work, an elongated box of dimensions 
D/ac — 16.7 and L/ac = 8.3 is used. 



Close to the critical point the simulation moves back 
and forth easily between the gas and liquid phases, while 
further away from the critical point, i.e. at higher poly- 
mer fugacity, the free energy barrier between the two 
phases increases. Hence transitions from one to the other 
phase become less likely, and the simulation spends most 
of the time in only one of the two phases. A crucial 
ingredient in our simulation is therefore the use of a 
biased sampling technique. We employ successive um- 
brella sampling, as was recently developed by Virnau and 
Miiller 29] , to enable accurate sampling in regions of rjc 
where P{ric), due to the free energy barrier separating 
the phases, is very small. In this approach, states (or 
windows) are sampled successively. In the first window, 
the number of colloids is allowed to vary between and 
1, in the second window between 1 and 2, and so forth. 
The number of polymers is allowed to fluctuate freely in 
each window. Our sampling scheme is thus strictly one- 
dimensional: the bias is put on the colloid density only. 
Note that this becomes problematic for very large sys- 
tems because of droplet formation. An additional free 
energy barrier, which grows with the size of the simu- 
lation box, must be crossed before the transition from 
the droplet state, to the slab geometry occurs (which is 
required if the interfacial tension is to be determined). 
As pointed out in Refs. fsil Is^ this additional barrier 
is not one in colloid density, but rather in the energy- 
like parameter, which for our system would be the poly- 
mer density. Hence, for very large systems, a naive one- 
dimensional biasing scheme such as described above, is 
prone to severe sampling difficulties (more appropriate in 
this case would be a two-dimensional sampling scheme in 
both the colloid and the polymer density). For the sys- 
tem size used by us, however, no problems in obtaining 
the slab geometry were encountered. 

In this work, we simulate up to a colloid packing frac- 
tion of rjc ~ 0.45, corresponding to a maximum of around 
1000 colloidal particles. The maximum number of poly- 
mers is obtained at low 77c. While the precise number 
depends on and /3e, a value of 3000 is typical. In 
each window, O{10^) grand canonical cluster moves are 
attempted, of which 0(10^) are accepted at low rjc, and 
O{10'^) at high rjc (the grand canonical cluster move thus 
becomes less efficient with increasing colloid packing frac- 
tion). The typical CPU time investment to obtain a sin- 
gle distribution P{ric) is 48 hours on a moderate com- 
puter. 



IV. RESULTS 

In Fig. n we plot results for the bulk fluid- fluid demix- 
ing binodal as obtained from theory and simulation in 
the (r^c, ??p)-plane, i.e. the "system representation". For 
the AOV case (recovered for /3e = 0), the DFT predicts 
the same bulk free energy for fluid states, and hence the 
same demixing binodal, as free volume theory [Tlj . The 
result is known to compare overall reasonably well with 
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simulation results for a variety of size ratios q, but to de- 
viate close to the critical point 0, 0, [s, 9] . For increasing 
strength of the polymer-polymer repulsion, the theoret- 
ical critical point shifts toward higher values of 77c, and 
very slightly to lower values of rjp. The accompanying 
shift of the binodal leads to a growth of the one-fluid 
region in the phase diagram, hence polymer-polymer re- 
pulsion tends to promote mixing. The simulation results 
indicate the same trend, but show a quantitatively larger 
shift of the binodal toward higher values of 77c upon in- 
creasing /3e. Note also the pronounced finite-size devia- 
tions of the simulation data in the vicinity of the critical 
point. As a result, the slight decrease in the critical value 
of rip with increasing f3e, as predicted by DFT, cannot be 
identified in the simulation data. To better access the 
critical region, a finite size scaling analysis [sj Is^. Issf 
would be required. While such an investigation has been 
carried out for the standard AOV model [sj, it would 
require extensive additional simulations for the extended 
model, which is beyond the scope of the present work. 

When plotted in the (r/c, ?7p)-plane or "reservoir repre- 
sentation" , the theoretical results display a similar shift 
of the binodal toward larger values of r]c, and consid- 
erable broadening of the coexistence region, see Fig. |2] 
This implies that at a given value of r]^ above the critical 
point, increasing the polymer-polymer repulsion leads to 
a larger density difference between the coexisting phases. 
The broadening of the coexistence region is strikingly 
confirmed by the simulation data. We again observe an 
overall shift of the binodals toward higher values of rjc- 
In accordance with findings for the standard AOV model 
Hills','?!, the critical value of r/^ obtained by DFT un- 
derestimates the simulation value. This is related to the 
universality class of the AOV model, which is that of the 
3D Ising model |3l)|, and hence, close to criticality, devi- 
ations from a mean-field treatment like the present DFT 
are to be expected. 

In Fig. 13 we show results for the gas-liquid interfa- 
cial tension, 7, plotted as a function of the difference 
between the colloid packing fraction in the coexisting 
phases, Aijc = Vc ~ Vi- Both DFT and simulation indi- 
cate that the interfacial tension decreases markedly upon 
increasing 0e, in accordance with square-gradient treat- 
ments |3^ I33 • Note that this change is not due to in- 
terfacial contributions alone, but also to the pronounced 
changes in the location of the bulk demixing binodal dis- 
cussed above. The decrease might seem reasonable based 
on the profound increase in the density of the coexisting 
liquid upon increasing /3e, see again the phase diagram 
in reservoir representation (Fig.|2J). As is apparent from 
Fig. 12 a given value of Aijc corresponds to a reduction 
of r^p at coexistence upon increasing (3e. We believe that 
this decrease induces the observed decrease in 7. 

Alternatively, comparing the interfacial tension as a 
function of rj^ thus reveals an increase with increasing 
/3e, see the upper panel of Fig. 0] which shows the DFT 
result. The lower panel of Fig. ^ shows the correspond- 
ing simulation data. The increase in the interfacial ten- 



sion with increasing /3e is confirmed up to /3e = 0.25, 
but not for /3e — 0.5. The simulation data for the lat- 
ter case, however, should be treated with some care, 
as in grand canonical simulations, the quantity rj^ does 
not play the role of a direct control parameter, which 
rather the polymer fugacity, Zp, does. Consequently, con- 
version to the reservoir representation requires an addi- 
tional step, namely the determination of rj^ as a func- 
tion of Zp. With the exception of /3e — 0, in which case 
T^p — TTCTpZp/G holds trivially, this conversion introduces 
an additional statistical uncertainty. Fig. [21 shows that, 
for /3e = 0.5, the binodal has become very flat, so even 
small uncertainties in r]p imply rather large uncertain- 
ties in Arjc- Hence, for very flat binodals, the reservoir 
representation is not convenient in simulations (more re- 
liable in this case is the system representation, see Fig. O 
which, as a benefit, is experimentally more relevant). A 
second point is that the Monte Carlo cluster move be- 
comes less efficient with increasing rjc and /3e (see Ref. 0) 
and this will also adversely affect the data (especially for 
f3e = 0.5, since then both /3e and yyj. are substantial). 
Therefore, we conclude that the shift of the simulation 
data for (3e — 0.5 in Fig. ^ most likely reflects a simula- 
tion artifact. 



V. CONCLUSIONS 

In conclusion, we have investigated the effect of 
polymer-polymer repulsion on the fluid-fluid demixing 
phase behavior and on the (colloidal) liquid-gas inter- 
face of a model colloid-polymer mixture. We have used a 
simplistic pair potential between polymers, given by a re- 
pulsive stepfunction, to extend the standard AOV model 
to cases of interacting polymers. Grand canonical Monte 
Carlo simulations of the full mixture demonstrate the rea- 
sonable accuracy of the theory, with a tendency to quan- 
titatively underestimate the shifts in the bulk fluid-fluid 
demixing binodal, and the gas-liquid interfacial tension, 
upon increasing strength of the polymer-polymer repul- 
sion. The present study demonstrates the usefulness of 
the model as such, as it clearly displays previously found 
features due to polymer non-ideality. This offers ways to 
study further interesting inhomogeneous situations, like 
the wetting properties at substrates. Such investigations 
could test the robustness of the results obtained for the 
surface pha se behavior at a hard wall using ideal poly- 
mers Hill in El m 
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FIG. 1: Bulk fluid-fluid demixing phase diagram of the ex- 
tended AOV model as a function of colloid packing frac- 
tion, ric and polymer packing fraction, rip, for size ratio 
q = 0.8 and increasing strength of the polymer-polymer repul- 
sion f3e — 0,0.1,0.25,0.5 as indicated, a) The binodal (lines) 
and critical point (symbols) as obtained from DFT. The case 
/3e = corresponds to the result from free volume theory for 
the AOV model, b) The binodaJ as obtained from simulations 
(symbols indicate data points; lines are guides to the eye). 
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FIG. 2: The analogue of Fig. d but as a function of colloid 
packing fraction, r^c, and polymer reservoir packing fraction 
rjp in a reservoir of interacting polymers. 
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FIG. 3: Dimensionless interfacial tension P-ycr^, of the free 
colloidal gas-liquid interface, as a function of the difference in 
colloid packing fractions of the coexisting states, Ar/c = r}]. — 
rjf. The size ratio is q = 0.8, and the strength of the polymer- 
polymer repulsion equals /3e — 0, 0.1, 0.25, 0.5 as indicated, a) 
Results from DFT; b) results from simulation. 
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FIG. 4: The analogue of Fig. H but as a function of the 
polymer reservoir packing fraction rip. a) Results from DFT; 
b) results from simulation. 



